Load in necessary packages
library(tidyverse)
## ── Attaching packages ──────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0 ✓ purrr 0.3.4
## ✓ tibble 3.0.1 ✓ dplyr 0.8.5
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ─────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(sf)
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(mapview)
library(tigris)
## To enable
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
##
## Attaching package: 'tigris'
## The following object is masked from 'package:graphics':
##
## plot
library(censusapi)
##
## Attaching package: 'censusapi'
## The following object is masked from 'package:methods':
##
## getFunction
library(leaflet)
library(lehdr)
library(fs)
options(
tigris_class = "sf",
tigris_use_cache = TRUE
)
Sys.setenv(CENSUS_KEY="10dcd73d7c043e91bac9fb8d3989cbff54b08790")
#Load in PUMAS
scc_pumas <- pumas("CA", cb=F, progress_bar=F) %>%
filter(grepl("Santa Clara County", NAMELSAD10)) #06 = CA, 085 = Santa Clara
scc_pumas_list <- scc_pumas %>%
pull(PUMACE10)
#Load in blockgroups
scc_blockgroups <- block_groups("CA","Santa Clara", cb=F, progress_bar=F)
mapview(scc_blockgroups)
scc_bg_list <- scc_blockgroups %>%
pull(GEOID)
#Get water area
water <-
"Santa Clara County" %>%
map_dfr(function(county){
area_water("CA", county, progress_bar = FALSE) %>% as.data.frame()
}) %>%
st_as_sf() %>%
st_set_crs(st_crs(counties("CA", cb=F, progress_bar = FALSE)
)) %>%
st_union() %>%
as_Spatial() %>%
sp::disaggregate() %>%
st_as_sf() %>%
st_set_crs(st_crs(counties("CA", cb=F, progress_bar = FALSE)
)) %>%
mutate(area = st_area(.) %>% as.numeric()) %>%
filter(area == max(area)) %>%
dplyr::select(-area)
## Warning in bind_rows_(x, .id): Vectorizing 'sfc_POLYGON' elements may not
## preserve their attributes
scc_pumas_fixed <- scc_pumas %>%
st_difference(water)
## although coordinates are longitude/latitude, st_difference assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
scc_blockgroups_fixed <- scc_blockgroups %>%
st_difference(water)
## although coordinates are longitude/latitude, st_difference assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
#Load in Santa Clara County Puma data
scc_puma_data <- readRDS(file = "/Users/spencerrobinson/Documents/covid19/Spencer_Robinson/bay_puma_soc.rds") %>% #PUMA = Public Use Microdata Area
filter(PUMA %in% scc_pumas_list) %>%
select("PUMA","SOCP","PWGTP")
scc_puma_data_weighted <- data.frame(PUMA = rep(scc_puma_data$PUMA, times = scc_puma_data$PWGTP), SOCP = rep(scc_puma_data$SOCP, times = scc_puma_data$PWGTP))
#Load in SOC data
soc_frac <- readRDS(file = "/Users/spencerrobinson/Documents/covid19/Spencer_Robinson/soc_fraction")
#Join SOC data to bay_puma_data for Santa Clara County only
scc_puma_data_soc <- scc_puma_data_weighted %>%
left_join(soc_frac, by = c("SOCP" = "SOC Occupation"))
## Warning: Column `SOCP`/`SOC Occupation` joining factor and character vector,
## coercing into character vector
# head(soc_essential_weighted_frac)
scc_puma_laborforce <- scc_puma_data_soc %>%
filter(!is.na(SOCP)) %>% #Remove those not in work force
group_by(PUMA) %>%
summarize(fracEssential = round(mean(fracEssential)*100, digits =1),population = n())
#
# head(scc_puma_data_soc)
#Now we'll visualize by PUMA
blue_pal <- colorNumeric(
palette = "Blues",
domain =
c(scc_puma_laborforce %>%
pull(fracEssential) %>%
unique())
)
puma_laborforce <- leaflet(data = scc_puma_laborforce) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
data =
scc_puma_laborforce %>%
left_join(scc_pumas_fixed, by = c("PUMA" = "PUMACE10")) %>%
st_as_sf() %>% st_set_crs(4326),
fillColor = ~blue_pal(fracEssential),
color = "white",
weight = 0.25,
fillOpacity = 0.5,
label = ~paste0(fracEssential,"%"),
highlightOptions = highlightOptions(
weight = 2,
opacity = 1
)
) %>%
addLegend(pal = blue_pal, values = ~fracEssential, opacity = 1, title = "% Essential Workers" )
## Warning: Column `PUMA`/`PUMACE10` joining factor and character vector, coercing
## into character vector
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
puma_laborforce
scc_pumas_w_area <- scc_pumas_fixed %>%
mutate(puma_area = st_area(.))
#Overlap PUMA1 and Blockgroup
mapview(scc_pumas_fixed, color = "red")+
mapview(scc_blockgroups_fixed, color = "blue")
#Try just with PUMA 1
bg_puma1 <- scc_pumas_w_area %>%
filter(PUMACE10 == "08501") %>%
st_intersection(scc_blockgroups_fixed) %>%
mutate(bg_area = st_area(.)) %>%
mutate(percent_area_puma = bg_area/puma_area)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
bg_puma1 %>% mapview()
bg_puma1 %>%
summarize(total = sum(percent_area_puma)) # Equals 1
## Simple feature collection with 1 feature and 1 field
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -122.2027 ymin: 37.20049 xmax: -122.0286 ymax: 37.46557
## CRS: 4269
## total geometry
## 1 0.9999996 [1] MULTIPOLYGON (((-122.044 37...
#PUMA intersect all SOCs
bg_puma <- scc_pumas_w_area %>%
st_intersection(scc_blockgroups_fixed) %>%
mutate(bg_area = st_area(.)) %>%
mutate(percent_area_puma = bg_area/puma_area)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
bg_puma <- bg_puma %>%
mutate(`bg_area` = as.double(bg_area)) %>%
filter(bg_area > 0)
#
# head(bg_puma)
bg_puma_laborforce <- bg_puma %>%
left_join(scc_puma_laborforce , by = c("PUMACE10" = "PUMA")) %>%
as_tibble() %>%
select(PUMACE10, GEOID,fracEssential) %>%
mutate(fracEssential = as.numeric(fracEssential))
## Warning: Column `PUMACE10`/`PUMA` joining character vector and factor, coercing
## into character vector
blue_pal <- colorNumeric(
palette = "Blues",
domain = c(bg_puma_laborforce %>%
pull(fracEssential)
)
)
bg_laborforce <- leaflet(data = bg_puma_laborforce) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
data = bg_puma_laborforce %>%
left_join(scc_blockgroups_fixed) %>%
st_as_sf() %>% st_set_crs(4326),
fillColor = ~blue_pal(fracEssential),
color = "white",
weight = 0.25,
fillOpacity = 0.5,
label = ~paste0(fracEssential,"%"),
highlightOptions = highlightOptions(
weight = 2,
opacity = 1
)
) %>%
addLegend(pal = blue_pal, values = ~fracEssential, opacity = 1, title = "% Essential of Labor Force" )
## Joining, by = "GEOID"
## Warning: st_crs<- : replacing crs does not reproject data; use st_transform for
## that
bg_laborforce
Social Distancing